Multivariate data analysis: a discriminating method of the origin of wines

source("https://raw.githubusercontent.com/MatteoFasulo/Rare-Earth/main/R/util/coreFunctions.R")

loadPackages(c('sn','tidyverse','psych','RColorBrewer','stargazer','mclust','ContaminatedMixt',
               'plotly','ggplot2','ggdendro','teigen','tclust','HDMD','caTools','clustvarsel',
               'vscc','sparcl','pgmm','caret','glmnet','MLmetrics','DT'))

load("Z:\\DesktopC\\LUMSA\\2\\Data Mining\\Finite Mixture\\FiniteMixtureL31.RData")
rm(CO2data)
rm(NOdata)
rm(tonedata)
type <- wine$Type
rm(wine)

1 Introduzione

Dati su 27 caratteristiche chimico/fisiche di tre diversi tipi di vino (Barolo, Grignolino, Barbera) dal Piemonte. Un set di dati con 178 osservazioni e 28 variabili (di cui la prima relativa alla tipologia di vino). Nell’ordine:

• Barolo
• Grignolino
• Barbera

data(wines)
wines

1.1 Annata del vino: un fattore da considerare?

E’ stato possibile attraverso la ricerca originaria risalire all’anno di osservazione di ciascun vino. Di seguito vengono riportate le osservazioni dei tre diversi tipi di vino durante gli anni:

year <- as.numeric(substr(rownames(wines), 6, 7))
table(wines$wine, year)

            year
             70 71 72 73 74 75 76 78 79
  Barolo      0 19  0 20 20  0  0  0  0
  Grignolino  9  9  7  9 16  9 12  0  0
  Barbera     0  0  0  0  9  0  5 29  5

#wines[,'wine'] <- type

Notiamo subito che il Barbera è distribuito principalmente negli ultimi anni (76,78,79) mentre il Barolo nel 71, 73 e 74. Per quanto riguarda la percentuale delle singole classi:

wines %>%
  count(wine = factor(wine)) %>%
  mutate(pct = prop.table(n)) %>% 
  ggplot(aes(x = wine, y = pct, fill = wine, label = scales::percent(pct))) + 
  geom_col(position = 'dodge') + 
  geom_text(position = position_dodge(width = .9),
            vjust = -0.5, 
            size = 3) + 
  scale_y_continuous(name = "Percentage")+
  scale_x_discrete(name = "Wine Name")+
  scale_fill_hue(l=40, c=35)+
  theme(legend.position = "none")

E’ chiaro che il Grignolino sia il più numeroso (39.9%) seguito dal Barolo (33.1%) e dal Barbera (27.0%).

Per rappresentare la dispersione dei dati abbiamo usato uno scatterplot leggermente differente dal solito. Sulla diagonale superiore si vede la distribuzione dei dati mentre sulla diagonale inferiore vi è la correlazione di Pearson tra le variabili.

# Correlation panel
my_cols <- c("#00AFBB", "#E7B800", "#FC4E07") 
panel.cor <- function(x, y){
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(0, 1, 0, 1))
    r <- round(cor(x, y), digits=2)
    txt <- paste0("R = ", r)
    text(0.5, 0.5, txt, cex = 1)
}
# Customize upper panel
upper.panel<-function(x, y){
  points(x,y, pch = 19, col = my_cols[wines$wine],cex=.5)
}
# Create the plots
pairs(wines[,2:8], 
      lower.panel = panel.cor,
      upper.panel = upper.panel)

Dalle prime 7 variabili, è possibile notare una correlazione di 0.69 tra acidity e malic e ovviamente una relazione inversamente proporzionale tra pH e acidity.

2 Statistiche descrittive

Per visualizzare le statistiche descrittive (media e deviazione standard) ci è sembrato opportuno dividerle in base alla classe di appartenenza:

printMeanAndSdByGroup <- function(variables,groupvariable)
  {
     variablenames <- c(names(groupvariable),names(as.data.frame(variables)))
     groupvariable <- groupvariable[,1]
     means <- aggregate(as.matrix(variables) ~ groupvariable, FUN = mean)
     names(means) <- variablenames
     print(paste("Means:"))
     print(means)
     sds <- aggregate(as.matrix(variables) ~ groupvariable, FUN = sd)
     names(sds) <- variablenames
     print(paste("Standard deviations:"))
     print(sds)
}
printMeanAndSdByGroup(wines[2:28],wines[1])

[1] "Means:"

[1] "Standard deviations:"

Alcune considerazioni:

• La media di sugar, potassium, magnesium, phosphate, chloride, flavanoids, proanthocyanins, colour nel Barolo è più alta.
• La media di acidity, tartaric, malic, uronic, alcal_ash nel Barbera è più alta.
• La deviazione standard di acidity è più alta nel Grignolino.

2.1 Varianza Within

Abbiamo calcolato la varianza within tra una feature e i tipi di vino:

calcWithinGroupsVariance(wines["flavanoids"],wines[1])

[1] 0.2747075

2.2 Varianza Between

Stesso discorso per la varianza between tra una feature e i vini:

calcBetweenGroupsVariance <- function(variable,groupvariable)
  {
     # find out how many values the group variable can take
     groupvariable2 <- as.factor(groupvariable[[1]])
     levels <- levels(groupvariable2)
     numlevels <- length(levels)
     # calculate the overall grand mean:
     grandmean <- sapply(variable,mean)
     # get the mean and standard deviation for each group:
     numtotal <- 0
     denomtotal <- 0
     for (i in 1:numlevels)
     {
        leveli <- levels[i]
        levelidata <- variable[groupvariable==leveli,]
        levelilength <- length(levelidata)
        # get the mean and standard deviation for group i:
        meani <- mean(levelidata)
        sdi <- sd(levelidata)
        numi <- levelilength * ((meani - grandmean)^2)
        denomi <- levelilength
        numtotal <- numtotal + numi
        denomtotal <- denomtotal + denomi
     }
     # calculate the between-groups variance
     Vb <- numtotal / (numlevels - 1)
     Vb <- Vb[[1]]
     return(Vb)
}
calcBetweenGroupsVariance(wines["flavanoids"],wines[1])

[1] 64.2612

2.3 Separazione

Per vedere quali variabili hanno la maggiore separazione, (rapporto tra Varianza Between e Varianza Within) abbiamo scritto una funzione apposita per calcolarne il valore per ogni feature.

3 Splitting in Train e Test

Per alcune analisi successive, abbiamo provato a cambiare il task della ricerca tentando di utilizzare un modello come classificatore e di misurarne le prestazioni di classificazione. A tale scopo, abbiamo suddiviso il dataset originario in due sottogruppi:

• train
• test

require(caTools)
sample = sample.split(wines[,1], SplitRatio = .50)

train = subset(wines, sample == TRUE)
trainTestNames <- train$wine
print(paste("Train Obs:",nrow(train)))

[1] "Train Obs: 90"

#train$wine <- as.numeric(train$wine)

test  = subset(wines, sample == FALSE)
wineTestNames <- test$wine
print(paste("Test Obs:",nrow(test)))

[1] "Test Obs: 88"

#test$wine <- as.numeric(test$wine)

4 Clustering

Per il Clustering abbiamo deciso di applicare:

• Un approccio Distance-Based:
  • Gerarchico:
    • Euclidean, Minkowski, Manhattan, Mahalanobis
  • Di partizionamento:
    • KMeans
    • PAM
• Un approccio Model-Based:
  • Gaussian Mixture (Mclust)
  • Contaminated Normal (CNmixt)
  • Multivariate t Distribution (teigen)
  • Parsimonious Gaussian Mixture Models (pgmm)

4.1 Distance-Based

dissMatrix    <- pairwise.mahalanobis(wines[,-1],
                                      grouping = c(1:nrow(wines)),
                                      cov = cov(wines[,-1]))$distance
dissMatrix <- sqrt(dissMatrix)
dissMatrix <- as.dist(dissMatrix)
combDist <- function(distance, methods, df, dt, dissMatrix) {
  c <- 0
  results <- list()
  for (i in 1:length(distance)){
    ifelse(distance[i] == "minkowski",
           dist <- dist(df, method = distance[i], p = 4),
           ifelse(distance[i] == "mahalanobis",
                  dist <- dissMatrix,
                  dist <- dist(df, method = distance[i])))
    for (j in 1:length(methods)){
      dendo = hclust(dist, method = methods[j])
      dendo.x = ggdendrogram(dendo, rotate = F, size = 2, leaf_labels = T, labels = F) + 
                               ggtitle(paste(distance[i],' ',methods[j],sep=''))
      for(elem in 2:4){
        cluster = cutree(dendo, k=elem)
        c <- c + 1
        results[[c]] <- list(distance = distance[i],
                             method = methods[j],
                             groups = elem,
                             table = table(dt,cluster),
                             dendo = dendo.x,
                             AdjustedRandIndex = adjustedRandIndex(dt,cluster),
                             cluster = cluster)
      }
    }
  }
  return(results)
}
results <- combDist(c("euclidean", "manhattan", "minkowski","mahalanobis"),
                    c("single", "complete", "average", "ward.D"), scale(wines[,-1]), wines[,1], dissMatrix)

optimal <- function(results){
  best_randIndex.eu = 0
  best_randIndex.ma = 0
  best_randIndex.mi = 0
  best_randIndex.maha = 0
  best_model.eu = integer()
  best_model.ma = integer()
  best_model.mi = integer()
  best_model.maha = integer()
  for (i in 1:length(results)){
    current_randIndex = results[[i]]$AdjustedRandIndex
    if (results[[i]]$distance == "euclidean"){
      if (current_randIndex > best_randIndex.eu) {
        best_randIndex.eu = current_randIndex
        best_model.eu = i
      }
    }
    else if (results[[i]]$distance == "manhattan"){
      if (current_randIndex > best_randIndex.ma) {
        best_randIndex.ma = current_randIndex
        best_model.ma = i
      }
    }
    else if (results[[i]]$distance == "minkowski"){
      if (current_randIndex > best_randIndex.mi) {
        best_randIndex.mi = current_randIndex
        best_model.mi = i
      }
    }
    else if (results[[i]]$distance == "mahalanobis"){
      if (current_randIndex > best_randIndex.maha) {
        best_randIndex.maha = current_randIndex
        best_model.maha = i
      }
    }
  }
  #print(list(euclidean = list(model.number = best_model.eu,
  #                            cluster = results[[best_model.eu]]$groups,
  #                            AdjustedRandIndex = best_randIndex.eu),
  #           manhattan = list(model.number = best_model.ma,
  #                            cluster = results[[best_model.ma]]$groups,
  #                            AdjustedRandIndex = best_randIndex.ma),
  #           minkowski = list(model.number = best_model.mi,
  #                            cluster = results[[best_model.mi]]$groups,
  #                            AdjustedRandIndex = best_randIndex.mi),
  #           mahalanobis=list(model.number = best_model.maha,
  #                            cluster = results[[best_model.maha]]$groups,
  #                            AdjustedRandIndex = best_randIndex.maha))
  #    )
  return(list(euclidean = results[[best_model.eu]],
              manhattan = results[[best_model.ma]],
              minkowski = results[[best_model.mi]],
              mahalanobis=results[[best_model.maha]]))
}
best_dist_model = optimal(results)

4.1.1 Dendrogrammi

4.1.1.1 Euclidean

print(paste("AdjustedRandIndex:",round(best_dist_model$euclidean$AdjustedRandIndex,3)))

[1] "AdjustedRandIndex: 0.769"

4.1.1.2 Manhattan

print(paste("AdjustedRandIndex:",round(best_dist_model$manhattan$AdjustedRandIndex,3)))

[1] "AdjustedRandIndex: 0.946"

4.1.1.3 Minkowski

ggplotly(best_dist_model$minkowski$dendo)

print(best_dist_model$minkowski$table)

            cluster
dt            1  2  3  4
  Barolo     55  2  2  0
  Grignolino  8  1 59  3
  Barbera     2  0  0 46

print(paste("AdjustedRandIndex:",round(best_dist_model$minkowski$AdjustedRandIndex,3)))

[1] "AdjustedRandIndex: 0.73"

4.1.1.4 Mahalanobis

print(paste("AdjustedRandIndex:",round(best_dist_model$mahalanobis$AdjustedRandIndex,4)))

[1] "AdjustedRandIndex: 0.27"

4.1.2 K-Means e PAM

4.1.2.1 KMeans 3

require(cluster)
k.means.3 <- kmeans(scale(wines[,-1]),centers=3,nstart = 50, iter.max = 100)
table(wines[,1], k.means.3$cluster)

            
              1  2  3
  Barolo      1  0 58
  Grignolino 68  1  2
  Barbera     0 48  0

print(paste("AdjustedRandIndex:",round(adjustedRandIndex(k.means.3$cluster, wines[,1]),3)))

[1] "AdjustedRandIndex: 0.93"

4.1.2.2 KMeans 4

require(cluster)
k.means.4 <- kmeans(scale(wines[,-1]),centers=4,nstart = 50, iter.max = 100)
table(wines[,1], k.means.4$cluster)

            
              1  2  3  4
  Barolo      0  0 59  0
  Grignolino 38  0  5 28
  Barbera     0 47  0  1

print(paste("AdjustedRandIndex:",round(adjustedRandIndex(k.means.4$cluster, wines[,1]),3)))

[1] "AdjustedRandIndex: 0.736"

4.1.2.3 PAM 3

require(cluster)
PAM.3 <- pam(wines[,-1], k=3,
    metric = "euclidean", 
    nstart = 50,
    stand = TRUE)
table(wines[,1], PAM.3$clustering)

            
              1  2  3
  Barolo     50  9  0
  Grignolino  3 67  1
  Barbera     1  1 46

print(paste("AdjustedRandIndex:",round(adjustedRandIndex(PAM.3$clustering, wines[,1]),3)))

[1] "AdjustedRandIndex: 0.754"

4.1.2.4 PAM 4

require(cluster)
PAM.4 <- pam(wines[,-1], k=4,
    metric = "euclidean", 
    nstart = 50,
    stand = TRUE)
table(wines[,1], PAM.4$clustering)

            
              1  2  3  4
  Barolo     28 27  4  0
  Grignolino  3  0 67  1
  Barbera     1  0  1 46

print(paste("AdjustedRandIndex:",round(adjustedRandIndex(PAM.4$clustering, wines[,1]),3)))

[1] "AdjustedRandIndex: 0.728"

5 Variable Selection

5.1 Principio filosofico

Novacula Occami: frustra fit per plura quod potest fieri per pauciora (Il rasoio di Occam: è futile fare con più mezzi ciò che si può fare con meno). Tale principio metodologico è ritenuto alla base del pensiero scientifico moderno.

5.1.1 Headlong

            
              1  2  3  4
  Barolo     47 10  2  0
  Grignolino  6  3 61  1
  Barbera     0  0  1 47
[1] "AdjustedRandIndex: 0.734"

5.1.2 Greedy

            
              1  2  3  4
  Barolo     47 10  2  0
  Grignolino  6  3 61  1
  Barbera     0  0  1 47
[1] "AdjustedRandIndex: 0.734"

5.1.3 VSCC

table(wines[,1], vscc.mclust$initialrun$classification) #Clustering results on full data set

            
              1  2  3
  Barolo     58  1  0
  Grignolino  1 70  0
  Barbera     0  2 46

print(paste("AdjustedRandIndex:",round(adjustedRandIndex(vscc.mclust$initialrun$classification, wines[,1]),3)))

[1] "AdjustedRandIndex: 0.931"

table(wines[,1], vscc.mclust$bestmodel$classification) #Clustering results on reduced data set

            
              1  2  3
  Barolo     59  0  0
  Grignolino  3 66  2
  Barbera     0  0 48

print(paste("AdjustedRandIndex:",round(adjustedRandIndex(vscc.mclust$bestmodel$classification, wines[,1]),3)))

[1] "AdjustedRandIndex: 0.913"

5.1.4 KMeansSparse 3

5.1.5 KMeansSparse 4

5.2 Model-Based

mixt.wines <- Mclust(wines[,-1],G=3:8)

cn.wines.mixt <- CNmixt(wines[,-1], G = 3, initialization = "mixt", seed = 1234, parallel = F, verbose = F)
cn.wines.kmeans <- CNmixt(wines[,-1], G = 3, initialization = "kmeans", seed = 1234, parallel = F, verbose = F)
cn.wines.rpost <- CNmixt(wines[,-1], G = 3, initialization = "random.post", seed = 1234, parallel = F, verbose = F)
cn.wines.rclass <- CNmixt(wines[,-1], G = 3, initialization = "random.clas", seed = 1234, parallel = F, verbose = F)


teigen.kmeans <- teigen(wines[,-1], Gs=3:4, init = 'kmeans', scale = T)
adjustedRandIndex(wines[,1],teigen.kmeans$classification)

teigen.classifier.kmeans <- teigen(train[,-1], Gs=3:4, init = 'kmeans', scale = T, known = train[,1])
predict(teigen.classifier.kmeans,test[,-1])
adjustedRandIndex(test[,1],predict(teigen.classifier.kmeans, test[,-1])$classification)

teigen.classifier.uniform <- teigen(train[,-1], Gs=3:4, init = 'uniform', scale = T, known = train[,1])
predict(teigen.classifier.uniform,test[,-1])
adjustedRandIndex(test[,1],predict(teigen.classifier.uniform, test[,-1])$classification)

teigen.hard <- teigen(wines[,-1], Gs=3:4, init = 'hard', scale = T)
adjustedRandIndex(wines[,1],teigen.hard$classification)
teigen.soft <- teigen(wines[,-1], Gs=3:4, init = 'soft', scale = T)
adjustedRandIndex(wines[,1],teigen.soft$classification)

6 Tecniche di regolarizzazione

Per la nostra analisi abbiamo voluto verificare l’efficienza di tre noti modelli di regolarizzazione attraverso il pacchetto caret:

• Ridge
• Lasso
• Elastic Net

lambda <- 10^seq(0, -2, length = 250)

6.1 Ridge

Il modello Ridge riduce i coefficienti, in modo che le variabili, con un contributo minore al risultato, abbiano i loro coefficienti vicini allo zero. Invece di forzarli a essere esattamente zero (come nel Lasso), li penalizziamo con un termine chiamato norma L2 costringendoli così a essere piccoli. In questo modo diminuiamo la complessità del modello senza eliminare nessuna variabile attraverso una costante chiamata lambda (\(\lambda\)) di penalizzazione: \[ L_{ridge}(\hat{\beta}) = \sum_{i = 1}^{n}{(y_i - x_i\hat{\beta})^2} + \lambda\sum_{k = 1}^{K}{\hat{\beta}_k^2} \]

# Build the model
set.seed(123)
ridge <- caret::train(
  x = train[,-1],
  y = factor(train[,1]),
  method = "glmnet",
  trControl = trainControl("cv", number = 10, classProbs = TRUE, summaryFunction = multiClassSummary),
  tuneGrid = expand.grid(alpha = 0, lambda=lambda),
  metric="Accuracy")
# Model coefficients
coef(ridge$finalModel, ridge$bestTune$lambda)
# Make predictions
predictions.ridge <- ridge %>% predict(test)
# Model prediction performance
tibble(
  trueValue = wineTestNames,
  predictedValue = predictions.ridge)

Il ridge è composto dalla somma dei residui quadrati più una penalità, definita dalla lettera Lambda, che è moltiplicata per la somma dei coefficienti quadrati \(\beta\). Quando \(\lambda = 0\), il termine di penalità non ha alcun effetto e il ridge produrrà i coefficienti minimi quadrati classici. Tuttavia, quando \(\lambda\) aumenta all’infinito, l’impatto della penalità aumenta e i coefficienti si avvicinano allo zero. Il ridge è particolarmente indicato quando si hanno molti dati multivariati con numero di feature maggiore del numero di osservazioni. Lo svantaggio, però, è che includerà tutti le feature nel modello finale, a differenza dei metodi di feature selection, che generalmente selezioneranno un insieme ridotto di variabili tra quelle disponibili.

caret::confusionMatrix(predictions.ridge, test$wine)

6.2 Lasso

Il Least Absolute Shrinkage and Selection Operator (LASSO) riduce i coefficienti verso lo zero penalizzando il modello con un termine di penalità chiamato norma L1, che è la somma dei coefficienti in valore assoluto: \[ L_{lasso}(\hat{\beta}) = \sum_{i = 1}^{n}{(y_i - x_i\hat{\beta})^2} + \lambda\sum_{k = 1}^{K}{|\hat{\beta}_k|} \]

# Build the model
set.seed(123)
lasso <- caret::train(
  x = train[,-1],
  y = factor(train[,1]),
  method = "glmnet",
  trControl = trainControl("cv", number = 10, classProbs = TRUE, summaryFunction = multiClassSummary),
  tuneGrid = expand.grid(alpha = 1, lambda=lambda),
  metric="Accuracy")
# Model coefficients
coef(lasso$finalModel, lasso$bestTune$lambda)
# Make predictions
predictions.lasso <- lasso %>% predict(test)
# Model prediction performance
tibble(
  trueValue = wineTestNames,
  predictedValue = predictions.lasso)

In questo caso la penalità ha l’effetto di forzare alcune delle stime dei coefficienti, con un contributo minore al modello, a essere esattamente uguale a zero. Il lasso, quindi, può anche essere visto come un’alternativa ai metodi di feature selection per eseguire la selezione delle variabili al fine di ridurre la complessità del modello.

Come nel ridge, è fondamentale selezionare un buon valore di \(\lambda\).

Quando lambda è piccolo, il risultato è molto vicino alla stima dei minimi quadrati. All’aumentare di lambda, si verifica una contrazione in modo da poter eliminare le variabili che sono a zero.

caret::confusionMatrix(predictions.lasso, test$wine)

6.3 Elastic Net

Elastic Net combina le proprietà di Ridge e Lasso penalizzando il modello usando sia la norma L2 che la norma L1: \[ L_{ElasticNet}(\hat{\beta}) = \frac{\sum_{i = 1}^{n}{(y_i - x_i\hat{\beta})^2}}{2n} + \lambda(\frac{1-\alpha}{2}\sum_{k = 1}^{K}{\hat{\beta}_k^2} + \alpha\sum_{k = 1}^{K}{|\hat{\beta}_k}|) \]

set.seed(123)
elastic <- caret::train(
  x = train[,-1],
  y = factor(train[,1]),
  method = "glmnet",
  trControl = trainControl("cv", number = 10, classProbs = TRUE, summaryFunction = multiClassSummary),
  metric="Accuracy")
# Model coefficients
coef(elastic$finalModel, elastic$bestTune$lambda)
# Make predictions
predictions.enet <- elastic %>% predict(test)
# Model prediction performance
tibble(
  trueValue = wineTestNames,
  predictedValue = predictions.enet)

Oltre a impostare e scegliere un valore lambda, l’elastic net ci consente anche di ottimizzare il parametro alfa dove \(\alpha = 0\) corrisponde a ridge e \(\alpha = 1\) al lasso.

Pertanto possiamo scegliere un valore \(\alpha\) compreso tra 0 e 1 per ottimizzare l’elastic net. Se tale valore è incluso in questo intervallo, si avrà una riduzione con alcuni portati a \(0\).

caret::confusionMatrix(predictions.enet, test$wine)

6.4 Summary delle prestazioni

models <- list(ridge = caret::confusionMatrix(predictions.ridge, test$wine), 
               lasso = caret::confusionMatrix(predictions.lasso,test$wine),
               elastic = caret::confusionMatrix(predictions.enet,test$wine))
#models %>% summary(metric = "Accuracy")

7 Parsimonious Gaussian Mixture Models

# Il rilassamento dato da (p-q)^2 > p+q è verificato per  1 <= q <= 20
pg.random.ccc <- pgmmEM(selectedWines[,-1], rG=3:4, rq=1:5, seed = 1234, relax = T)
adjustedRandIndex(wines[,1],pg.random.ccc$map)

pg.random <- pgmmEM(wines[,-1], rG=3:4, rq=20:25,icl=F,zstart = 1,cccStart = F,loop = 50, seed = 1234, relax = T)
adjustedRandIndex(wines[,1],pg.random.ccc$map)

pg2 <- pgmmEM(wines[,-1], rG=3:4, rq=20:25,icl=T,zstart = 1,cccStart = T,loop = 50, seed = 1234, relax = T)
adjustedRandIndex(wines[,1],pg2$map)

pg2 <- pgmmEM(wines[,-1], rG=3:4, rq=20:25,icl=T,zstart = 2, seed = 1234, relax = T)
adjustedRandIndex(wines[,1],pg2$map)

pg2 <- pgmmEM(wines[,-1], rG=3:4, rq=20:25,icl=T,zstart = 3, seed = 1234, relax = T, zlist=)
adjustedRandIndex(wines[,1],pg2$map)




agree(cn.wines.mixt, givgroup = wines[,1])
agree(cn.wines.kmeans, givgroup = wines[,1])
agree(cn.wines.rpost, givgroup = wines[,1])
agree(cn.wines.rclass, givgroup = wines[,1])

plot(cn.wines.mixt, contours = TRUE)
plot(cn.wines.kmeans, contours = TRUE)
plot(cn.wines.rpost, contours = TRUE)
plot(cn.wines.rclass, contours = TRUE)

if (class(models[[i]])==‘kmeans’){ ariTrain[i] <- adjustedRandIndex(train,models[[i]]\(cluster) ariTest[i] <- adjustedRandIndex(test,models[[i]]\)cluster) bic[i] <- NA icl[i] <- NA G[i] <- length(unique(models[[i]]$cluster))

summaryOfModels <- function(train, test, models){
  nModel <- length(models)
  ariTrain <- NULL
  ariTest <- NULL
  bic <- NULL
  icl <- NULL
  G <- NULL
  for(i in 1:nModel){
           if (class(models[[i]])=='Mclust'){
      ariTrain[i] <- adjustedRandIndex(train,models[[i]]$classification)
      ariTest[i] <- adjustedRandIndex(test,predict(models[[i]],test)$classification)
      bic[i] <- models[[i]]$bic
      icl[i] <- models[[i]]$icl
      G[i] <- models[[i]]$G
    } else if (class(models[[i]])=='ContaminatedMixt'){
      bestModel <- getBestModel(models[[i]])$models[[1]]
      ariTrain[i] <- adjustedRandIndex(train,bestModel$group)
      ariTest[i] <- adjustedRandIndex(test,predict(bestModel,))
      bic[i] <- whichBestModel(models[[i]],train)$BIC
      icl[i] <- whichBestModel(models[[i]],train)$ICL
      G[i] <- whichBestModel(models[[i]],train)$G
    } else if (class(models[[i]])=='teigen'){
      ariTrain[i] <- adjustedRandIndex(train,models[[i]]$iclresults$classification)
      ariTest[i] <- adjustedRandIndex(test,models[[i]]$iclresults$classification)
      bic[i] <- models[[i]]$bic
      icl[i] <- models[[i]]$iclresults$icl
      G[i] <- models[[i]]$G
    } else if (class(models[[i]])=='pgmm'){
      ari[i] <- adjustedRandIndex(df,models[[i]]$map)
      G[i] <- models[[i]]$g
      ifelse(is.null(models[[i]]$bic[1])==TRUE,
             bic[i] <- NA,
             bic[i] <- as.double(models[[i]]$bic[1]))
      ifelse(is.null(models[[i]]$icl[1])==TRUE,
             icl[i] <- NA,
             icl[i] <- as.double(models[[i]]$icl[1]))
    } else if (class(models[[i]])=='tkmeans'){
      ari[i] <- adjustedRandIndex(df,models[[i]]$cluster)
      G[i] <- models[[i]]$k
      bic[i] <- NA
      icl[i] <- NA
    }
  }
  outputDF <- data.frame('AdjustedRandIndex' = ari,
                         'BIC' = bic,
                         'ICL' = icl,
                         'G' = as.integer(G),
                         row.names = c('KMeans','Mclust','ContaminatedMixt','TEigen','PGMM','TClust'),
                         stringsAsFactors = F)
  return(outputDF)
}
test <- summaryOfModels(wines[,1],list(k.means,mixt.wines,mixt.cn.wines,teigen_wine,pippo,trimmedOne))
test